#include <stdio.h>
#include "optimization.h"

using namespace alglib;


//Task parameters
struct TaskParams {
	double leftValue; //Left boundary value
	double rigthValue; //Right boundary value
	double nuLeft; //Right point coordinate
	double nuRight; //Left point coordinate
	int NCells; //Number of cells in grid
	double* nu; //Grid points coordinates
	double dNu; //Grid step
};

//Right hand side
double R(double x) {
	return 0;
};

//Target function
void  FRes(const real_1d_array &fi, real_1d_array &F, void *ptr)
{	
    TaskParams* taskParams = (TaskParams*)ptr;
	//double F[101];
	for (int i = 0; i<=taskParams->NCells; i++) {
		//Boundary conditions
		if (i == 0) {
			F[i] = fi[i] - taskParams->leftValue;
			continue;
		}
		if (i == taskParams->NCells) {
			F[i] = fi[i] - taskParams->rigthValue;
			continue;
		}

		//Compute derivatives
		//Second
		double d2 = (fi[i-1] - 2*fi[i] + fi[i+1]) / (taskParams->dNu * taskParams->dNu);

		//Third
		double d3 = 0;		
		if (i == 1) {
			double f0 = -fi[i];
			double f1 = +3*fi[i+1];
			double f2 = -3*fi[i+2];
			double f3 = 1 * fi[i+3];
			d3 = f0 + f1 + f2 + f3;
		} else if (i == taskParams->NCells-1) {
			double f0 = fi[i];
			double f1 = -3*fi[i-1];
			double f2 = +3*fi[i-2];
			double f3 = -1 * fi[i-3];
			d3 = f0 + f1 + f2 + f3;
		} else {
			double f0 = (fi[i-2] + fi[i-1]) / 2.0;
			double f1 = (fi[i-1] + fi[i]) / 2.0;
			double f2 = (fi[i] + fi[i+1]) / 2.0;
			double f3 = (fi[i+1] + fi[i+2]) / 2.0;
			d3 = -f0 + 3*f1 - 3*f2 + f3;
		};

		//Equation to solve
		F[i] = d3 - fi[i]*d2 - R(taskParams->nu[i]);
	};
	return;
}

//Main program ))
int main_(int argc, char *argv[]) {
	//Task parameters
	TaskParams taskParams;
	taskParams.NCells = 100; //Number of cells
	taskParams.nuLeft = 0.0;
	taskParams.nuRight = 5.0;
	taskParams.leftValue = 0.0;
	taskParams.rigthValue = 0.0;
	taskParams.dNu = (taskParams.nuRight - taskParams.nuLeft) / taskParams.NCells;

	//Set computational grid
	taskParams.nu = new double[taskParams.NCells+1];
	double deltaNu = taskParams.dNu;
	for (int i = 0; i<=taskParams.NCells; i++) {
		taskParams.nu[i] = taskParams.nuLeft + deltaNu * i;
	};	

	//Set initial guess
	real_1d_array fi;
	fi.setlength(taskParams.NCells+1); 
	for (int i = 0; i<=taskParams.NCells; i++) {
		fi[i] = 1;
	};	

	//Set optimization parameters
    double epsg = 0.0000000001;
    double epsf = 0;
    double epsx = 0;
    ae_int_t maxits = 1000;
    minlmstate state;
    minlmreport rep;

	//Run optimization
    minlmcreatev(taskParams.NCells+1, fi, 0.0001, state);
    minlmsetcond(state, epsg, epsf, epsx, maxits);
    minlmoptimize(state, FRes, NULL, (void *)&taskParams);
    minlmresults(state, fi, rep);

	//Output result
    printf("%d\n", int(rep.terminationtype)); // EXPECTED: 4
    printf("%s\n", fi.tostring(20).c_str()); // 

	FILE* fout = fopen("res.dat", "w");
	for (int i = 0; i<=taskParams.NCells; i++) {
		fprintf(fout, "%lf %lf\n", taskParams.nu[i], fi[i]); // 
	};
	fclose(fout);
    return 0;
}